home *** CD-ROM | disk | FTP | other *** search
/ The PC-SIG Library 10 / The PC-Sig Library - Shareware for the IBM PC and Compatibles (PC-SIG)(Tenth Edition Disks 1-2804)(1991).iso / PC_SIGCD / 04 / 2 / DISK0425.ZIP / ERFSIMP.PAS < prev    next >
Pascal/Delphi Source File  |  1985-04-03  |  2KB  |  79 lines
  1. program erfsimp;        { -> 323 }
  2.  
  3. { integration by Simpson's method }
  4.  
  5. const    tol        = 1.0E-4;
  6.     pi        = 3.141593;
  7.  
  8. var    done        : boolean;
  9.     sum,upper,lower,
  10.     erf,twopi    : real;
  11.  
  12. external procedure cls;
  13.  
  14. function fx(x: real): real;
  15. begin
  16.   fx:=exp(-x*x)
  17. end;        { function fx }
  18.  
  19.  
  20. procedure simps(function f(x: real): real;
  21.         lower,upper,tol    : real;
  22.         var sum        : real);
  23.  
  24. { numerical integration by Simpson's rule }
  25. { function is f, limits are lower and upper }
  26. { with number of regions equal to pieces }
  27. { partition is delta_x, answer is sum }
  28.  
  29. var    i        : integer;
  30.     x,delta_x,even_sum,
  31.     odd_sum,end_sum,
  32.     sum1        : real;
  33.     pieces        : integer;
  34. begin
  35.   pieces:=2;
  36.   delta_x:=(upper-lower)/pieces;
  37.   odd_sum:=f(lower+delta_x);
  38.   even_sum:=0.0;
  39.   end_sum:=f(lower)+f(upper);
  40.   sum:=(end_sum+4.0*odd_sum)*delta_x/3.0;
  41.   writeln(pieces:5,sum);
  42.   repeat
  43.     pieces:=pieces*2;
  44.     sum1:=sum;
  45.     delta_x:=(upper-lower)/pieces;
  46.     even_sum:=even_sum+odd_sum;
  47.     odd_sum:=0.0;
  48.     for i:=1 to pieces div 2 do
  49.       begin
  50.     x:=lower+delta_x*(2.0*i-1.0);
  51.     odd_sum:=odd_sum+f(x)
  52.       end;
  53.     sum:=(end_sum+4.0*odd_sum+2.0*even_sum)*delta_x/3.0;
  54.   until abs(sum-sum1)<=abs(tol*sum1)
  55. end;    { simps }
  56.  
  57. begin        { main program }
  58.   cls;
  59.   done:=false;
  60.   twopi:=2.0/sqrt(pi);
  61.   lower:=0.0;
  62.  
  63.   repeat
  64.     writeln;
  65.     writeln('Erf? ');
  66.     readln(upper);
  67.     if upper<0.0 then done:=true
  68.     else if upper=0.0 then
  69.       writeln('Erf of 0.0 is 0.0')
  70.  
  71.     else    { upper>0 }
  72.       begin
  73.         simps(fx,lower,upper,tol,sum);
  74.         erf:=twopi*sum;
  75.         writeln('Erf of ',upper:7:2,', is ',erf:8:4)
  76.         end
  77.     until done
  78. end.
  79.